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Abstract 

We present a general method for predicting bond percolation thresholds and critical surfaces 
for a broad class of two-dimensional periodic lattices, reproducing many known exact results and 
providing excellent approximations for several unsolved lattices. For the checkerboard and inho- 
mogeneous bow-tie lattices, the method yields predictions that agree with numerical measurements 
to more than six figures, and are possibly exact. 



*Electronic address: scullard@uchicago.edu 
^Electronic address: rziff@umich.edu 



1 



Percolation deals with the formation of long-range connectivity in a random network, 
represented by a given graph. In bond percolation, we define every edge of a graph to be 
open with probability p and closed with probability 1 — p. The average size of resulting 
clusters of open bonds becomes greater with increasing p and at a special value, p c , called 
the critical threshold, an infinite cluster appears. The value of p c strongly depends on the 
lattice under consideration. 

Knowledge of the percolation thresholds is central to understanding and applying this 
process, and the problem of deriving or computing those thresholds has been the subject 
of intense study since the introduction of the model over 50 years ago . Even after so 
many years of research there still exists no general method for deriving critical probabilities of 
arbitrary graphs, although schemes of varying accuracy have been proposed 0, Q, Ql . All non- 
trivial exact solutions are confined to two dimensions and special lattice classes p, Q] ; these 
exact thresholds all appear as the roots of polynomials with integer coefficients. Although 
we focus mainly on the Archimedean lattices in this Letter (Fig. [I]), we present a method 
whereby a polynomial can, in principle, be associated with any two-dimensional periodic 
lattice and the prediction for the bond threshold found as the root in [0, 1]. Our results 
agree with all exact solutions, and for unsolved problems we find excellent approximations 
for every system considered. In two cases, the checkerboard and inhomogeneous bow-tie 
lattices, we find predictions that agree with our numerical tests to high precision and these 
formulas may be exact. 



All known exact thresholds follow from duality arguments. Using a general triangle- 
triangle transformation [6J, which is an extension of the well-known star-triangle transfor- 
mation 181, broad classes of exact thresholds and critical surfaces have been derived recently 

n 

[9|]. The approach works only for lattices in which the basic repeated cell can be contained 
in a triangle, and the triangles are arranged in a self-dual way, as shown for example in Fig. 
UK Referring to the generalized cell in Figure [2b, the critical threshold is found by using 



the condition 



P(A,B,C) = P(A,B,C) (1) 
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FIG. 1: The Archimedean lattices; a) square (4 4 ); b) honeycomb (6 3 ); c) triangular (3 6 ); d) kagome 
(6,3,6,3); e) (4,8 2 );f) (3,12 2 );g) (3 3 ,4 2 ); h) (3,4,6,4); i) (3 2 , 4, 3, 4); j) (3 4 ,6);k) (4,6,12) . 
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FIG. 2: a) A class of exactly solvable lattices; b) every shaded triangle of a) can represent any 
network of sites and bonds contained in the vertices (A, B, C). 

where P(A, B, C) is the probability that the vertices (A, B, C) are connected by open bonds 
and P(A, B, C) is the probability that none are connected. As a simple example, we take for 
Figure Wp the star with inhomogeneous probabilities assigned as in Figure [3J which results 
in the honeycomb lattice (Figure db). By allowing distinct probabilities for the bonds on 
the unit cell, we find a critical surface as opposed to a critical point, and application of ([I]) 
gives jg] 

H(p, r, s) = prs — rp — rs — ps + 1 = 0. (2) 
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FIG. 4: a) The martini-A lattice; b) The assignment of probabilities on the unit cell. 

We locate the critical probability for the homogeneous system by setting r = s = p: 

p 3 -3p 2 + l = (3) 

so p c = 1 — 2sin7r/18 « 0.652704 [8]. If we want to find the threshold of the "martini-A" 
lattice (Fig. 0^,), we use the A cell shown in Figure H)d, which leads to 

A(pi,p2,ri,r2,r 3 ) = 1 - pir 2 - p 2 n - p\p 2 r 3 
- pinr 3 - p 2 r 2 r 3 + pip 2 r x r 3 
+ Pip 2 r 2 r 3 + p 1 r 1 r 2 r 3 

+ p 2 r 1 r 2 r 3 - Pip 2 r l r 2 r 3 = 0. (4) 



Setting all probabilities equal leads to the homogeneous polynomial and the bond threshold 
p c = 0.625457... [6]. Problems that can be solved with equation ([T]) have critical surface 
formulas of the same basic form as these two examples. For a lattice of the type in Figure [2] 
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FIG. 5: a) The assignment of probabilities on the (4, 8 2 ) unit cell; b) A "re-partitioning" of the 
graph into different unit cells. 

in which each of the n bonds of the unit cell is assigned a different probability p it this form 
is 

f(pi,P2,-,Pn) = 

where the function / is at most first-order in each argument. We refer to this latter property 
of / as "linearity". Since equation (pQ) compares events that happen only on a unit cell, it 
is clear that any / found through ([T]) will have this property. A consequence of linearity is 
that the critical threshold for a homogeneous system is the root of a polynomial of order at 
most n. 

To construct the critical surface for a general periodic lattice, our basic assumption is 
that linearity holds in general, even for lattices that cannot be solved with ([1]). In order 
to find / for a given lattice, we simply write the most general function that is first order 
in all its arguments, and then impose symmetries and agreement with special cases until 
the unknown coefficients are constrained. As an illustrative example, we choose the (4, 8 2 ) 
lattice of Figure [1^ with the probabilities assigned on the unit cell as in Figure [5^. The 
critical function can be written 



F(p, r,s,t,u,v) = aijkimnP l r J s k t 



l u m v n 



(5) 



i=0 n=0 



where there are 64 o's to be determined. In general, with a unit cell of n bonds there are 
2 n coefficients. Since can be multiplied by an arbitrary factor without changing the 
threshold, we are free to choose one of the o's, and we set the constant term to +1 as a 
matter of convention. To determine the rest, we impose symmetries and consider special 
cases. Since there are so many coefficients, we must use a computer algebra program to 
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impose the various constraints. We start with the repartitioning symmetry shown in Figure 
[5b. That is, there is more than one way to divide the (4, 8 2 ) lattice into unit cells and 
different choices should not alter the critical function, so 

F(p,r,s,t,u,v) = F(r,t,p,s,v,u). (6) 

This constraint requires many of the a's to be equal and reduces the number to be determined 
to 21. Next, we impose mirror symmetry; reflecting the lattice about a vertical line through 
the center of the unit cell does not change the threshold. This means 

F(p, r, s, t, u, v) = F(p } s, r, t, v, u) 

which reduces the number of coefficients by 3, to 18. Using a different repartitioning besides 
(EJ) and reflecting about a horizontal line do not give any additional constraints. We now 
turn to special cases. If we remove the t bond by setting its probability to we are left with 
a honeycomb lattice with two doubled bonds. Therefore, we require 

F(p, r, s, 0, u, v) = H{p, ru, sv) 
which leaves only three undetermined a's. If we set p = 1 we get the martini-A lattice, so 

F(l, r, s, t, u, v) = A(s, r, t, u, v) 

which fixes the remaining coefficients. We finally arrive at the critical surface, 

F(p, r, s, t,u,v) = 1 — pru — stu — psv — rtv 

— rsuv — ptuv + prsuv + ptruv 
+ pstuv + rstuv + prstu 

+ prstv — 2prstuv = 0. (7) 

The polynomial for the homogeneous case is 

F(p, p, p, p,p,p) = l- Ap 3 - 2p 4 + Qp 5 - 2p 6 = (8) 

with the solution p c = 0.676835... on [0, 1]. Parviainen's numerical estimate for this thresh- 
old is p c = 0.67680232(63), where the number in parentheses is the standard error in the last 
digit (s). Although our solution is ruled out, it is accurate to four significant figures. Addi- 
tionally, we have an analytic expression for the full critical surface, something that would 
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TABLE I: Comparison of bond percolation estimates with numerical results and the Riordan and 
Walters confidence intervals; a: b: [ill ], c: [lo| . d: 1 — 3p 2 — 6p 3 + 12p 4 — 6p 5 + p 6 = 0, e: 



1 - 3p 4 - 6p b + 3p 6 + 15p Y - 15/ + 4/ = 0, f: Eq. 
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;: l - 4p 2 - 12/ + 104/ - 193/ + U6p< 



45/ + 2p 10 = 0, h: 1 - 6/ - 12/ - 6/ + 69/ + 60/ - 363/ + 448/ - 252/° + 66p 
i: 1 - 4/ - 12p 3 - 2/ + 106p 5 - 186/ + 132/ - 36/ - 2p 9 + 2p 10 = . 



11 6p 12 = 0, 



be practically impossible to obtain numerically. While there are alternative paths to derive 
([7|) that are more efficient, the steps presented here demonstrate the general approach used 
for more complicated lattices. In fact, we have employed it to obtain critical surfaces and 
polynomials for all but two of the Archimedean lattices. These results are shown in Table [I] 

n 

along with the rigorous confidence intervals of Riordan and Walters [10] and the numerical 
values of Parviainen 



conjecture of Wu 



111 ]. Note that our estimate for the kagome lattice matches the 1979 
12], who took a similar approach in an attempt to solve that problem. 
Although all of our thresholds for the Archimedean lattices agree to at least four sig- 
nificant figures with Parviainen's results, they are ruled out numerically, implying that the 
linearity hypothesis does not hold exactly in these cases. However, it is possible to get 
successively better estimates with this method by increasing the number of distinct prob- 
abilities. Although there is a minimum number that one can choose, corresponding to the 
size of the unit cell, if we extend the inhomogeneity to a neighboring cell we get a refinement 
of the approximation. For the (3 3 , 4 2 ) lattice for example, the simplest unit cell contains 
five bonds, and using five probabilities leads to a polynomial (1 — 2p — 2p 2 + 3p 3 — /) with 
root in [0, 1], p c = 0.419308... . However, the value in Tabled] is based on covering two unit 
cells by using ten probabilities, giving a tenth-order polynomial. This estimate agrees much 



better with the numerical result and falls within the confidence interval of Ref. 



10(| . In gen- 
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eral, however, it becomes increasingly difficult to constrain systems with greater numbers of 
distinct probabilities, since the set of coefficients quickly becomes large. As such, all other 
results in Tabled] are based on a single unit cell. In fact, we have yet to find expressions for 
even the single-cell critical surfaces of the (3 4 , 6) and (4, 6, 12) lattices because of the large 
numbers of bonds involved. 

Our method also gives some solutions that agree with simulations to very high precision, 
and may be exact. If we set u — v — 1 in (jTJ) we get the square lattice (see Fig. [6^) but 
with four different probabilities. Equation ([7]) predicts 

C (p, r,s,t) = 1 — pr — ps — rs — pt — rt — st 

+ prs + prt + pst + rst = 0. (9) 

Wu [14] considered this situation, which he called the "checkerboard" , in the context of the 
q— state Potts model, for which percolation is the limit as q — > 1. Equation (Q is equivalent 
to his formula in that limit. Although Wu's conjecture was shown to be incorrect for q = 3 

n 

[15] . it holds in the q = 2 (Ising) case, and its validity for q = 1 (percolation) evidently 
has not been tested. We have carried out high-precision numerical simulations, using the 

n 

gradient percolation method described in [16], which are consistent with ([§]) being exact. To 
test predictions of the critical surface, we fix some probabilities and find the critical value 
for a single remaining parameter. For example, if we set p = 73/90, and assume t = s = r, 
then ([9]) predicts r c = 0.4. Numerically, we find r c = 0.40000004(10) simulating a total of 
10 14 frontier bonds in lattices of various gradients. We also found similar agreement at other 
parameter values. Note that (IHl) cannot be derived from condition ([T]). 



If we set v = 1 in ([7]), we get a prediction for the dual of the bow-tie lattice [l7|] (Figure [7J. 
The bond thresholds of dual pairs of lattices are related |8| by p c (L) = 1 — p c (L d ). This can 
be easily generalized for critical surfaces. Specifically, we find the 5-bond version (Figure 
[6b) of the bow-tie critical surface as 



F(l - p, 1 - r, 1 - s, 1 - t, 1 - u, 1) 



0. 
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FIG. 6: a) The inhomogeneous square lattice with the checkerboard assignment of probabilities; 
b) the assignment of probabilities on the 5-bond inhomogeneous bow-tie lattice (see Figure [7]). 




FIG. 7: The duality transformation for the bow-tie lattice (solid lines). The dotted graph is the 
dual. 

We thus find, for the lattice of Figure [6b, 

B(p, r,s,t,u) = 1 — u — pr — ps — rs — pt 
— rt — st + prs + prt + pst 
+ rst + pru + stu — prstu = 

which reduces to (j5J) in the limit u = 0. Setting the probabilities equal yields 

B(p, p, p, p, p) = 1 — p — 6p 2 + 6p 3 — p 5 = 



(10) 



with solution p c = 0.404518..., in agreement with Wierman's exact result, [17j]. We studied 
the situation in Figure [6)d numerically and once again find agreement to high precision. For 
example, if we set u = r = s = p and t = 1/2, then (fit)]) gives p c = (3 — v / 5)/2 = 0.38196601. 
Numerically, we find p c = 0.3819654(5), which is well within two standard deviations of the 
prediction. Like (jT0"|) cannot be derived using ([T]). 

Thus, in conclusion, we have a method that gives very accurate approximations to the 
bond critical surfaces for many lattices, and in two cases - the checkerboard and inho- 



mogeneous bow-tie lattices - formulas that are consistent with precise numerical results. 
Although the procedure presented here would work in principle for site percolation, there 
are not many exact site critical thresholds to which the Archimedean lattices can reduce in 
special cases. This problem also limits the extension to higher dimensions, where there are 
no exact solutions of any kind. 

RMZ acknowledges support of the National Science Foundation under grant DMS- 
0553487. 



[1] S. R. Broadbent and J. M. Hammersley, Proc. Cambridge Phil. Soc. 53, 629 (1957). 

[2] P. J. Flory, J. Amer. Chem. Soc. 63, 3083 (1941). 

[3] S. Galam and A. Mauger, Phys. Rev. E 53, 2177 (1996). 

[4] J. C. Wierman, D. P. Naor, and R. Cheng, Phys. Rev. E 72, 066116 (2005). 

[5] R. A. Neher, K. Mecke, and H. Wagner, J. Stat. Mech. p. P01011 (2008). 

[6] R. M. Ziff, Phys. Rev. E 73, 016134 (2006). 

[7] L. Chayes and H. K. Lei, J. Stat. Phys. 122, 647 (2006). 

[8] M. F. Sykes and J. W. Essam, J. Math. Phys. 5, 1117 (1964). 

[9] R. M. Ziff and C. R. Scullard, J. Phys. A 39, 15083 (2006). 

[10] O. Riordan and M. Walters, Phys. Rev. E 76, 011110 (2007). 

[11] R. Parviainen, J. Phys. A 40, 9253 (2007). 

[12] F. Y. Wu, J. Phys. C 12 (1979). 

[13] R. M. Ziff and P. N. Suding, J. Phys. A 15, 5351 (1997). 

[14] F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982). 

[15] I. G. Enting, J. Phys. A 20, L917 (1987). 

[16] R. M. Ziff and B. Sapoval, J. Phys. A 19, L1169 (1986). 

[17] J. C. Wierman, J. Phys. A 17, 1525 (1984). 



10 



